function WhiskerFit(varargin)
  
P = parsePairs(varargin);
checkField(P,'Frame');
checkField(P,'PStart',[45,180,40,100,20,40])
checkField(P,'FIG',1);
global cFrameFitting; 
cFrameFitting = double(P.Frame);

% PREPROCESSING WITH A WAVELET TYPE OF FILTER
% F = zeros(11); S = 1; 
% for iD = 0:35; 
%   for iX = -5:5; 
%     for iY=-5:5 
%       R = 2*pi/360*(iD*10); 
%       L = sqrt(iX^2+iY^2); 
%       if L<=10; 
%         V = [cos(R),sin(R)]'; 
%         O = [iX,iY]*V./(V'*V); 
%         D = sqrt(L^2-O^2); 
%         F(iX+11,iY+11,iD+1) = exp(-D^2/S^2)*cos(D*2); 
%       end; 
%     end; 
%   end; 
% end;
% for i=1:36  R(:,:,i) = conv2(E,F(:,:,i)); end
% RM = max(R,[], 3);

% FIT WHISKER
Opt = optimset; Opt.DiffMinChange = 1;
PFinal = lsqnonlin(@BezierFit,P.PStart,[],[],Opt);

Bezier = @(T,P0,P1,P2)bsxfun(@times,(1-T),bsxfun(@times,(1-T),P0)+bsxfun(@times,T,P1))...
  + bsxfun(@times,T,bsxfun(@times,(1-T),P1)+bsxfun(@times,T,P2));

% PLOT RESULTS
figure(P.FIG); clf; colormap(gray(256));
imagesc(cFrameFitting); set(gca,'yDir','normal'); hold on; caxis([-10,10]); 
T = [0:0.01:1];
BCurveInit= Bezier(T,P.PStart(1:2)',P.PStart(3:4)',P.PStart(5:6)');
plot(BCurveInit(1,:),BCurveInit(2,:),'r--','LineWidth',2);
BCurve= Bezier(T,PFinal(1:2)',PFinal(3:4)',PFinal(5:6)');
plot(BCurve(1,:),BCurve(2,:),'r-','LineWidth',2);

function E = BezierFit(P)

global cFrameFitting
P0 = P(1:2)';
P1 = P(3:4)';
P2 = P(5:6)';
T = [0:0.01:1];
Curve = bsxfun(@times,(1-T),bsxfun(@times,(1-T),P0) +bsxfun(@times,T,P1)) ...
  + bsxfun(@times,T,bsxfun(@times,(1-T),P1)+bsxfun(@times,T,P2));
B = zeros(size(cFrameFitting));
for i=1:length(T) 
  cM = Curve(:,i);
  B(round(cM(2)),round(cM(1))) = 20;
end
B = B(1:size(cFrameFitting,1),1:size(cFrameFitting,2));
B = relaxc2(B,1);
E = cFrameFitting - B;
E = E(:);










